#pragma once
#include <assert.h>
#include <stdio.h>
#include <string.h>
#include <time.h>

#include <algorithm>
#include <complex>
#include <cstdio>
#include <list>
#include <numeric>
#include <set>
#include <unordered_map>
#include <unordered_set>
#include <vector>

class Util
{
private:
    static constexpr double ZERO_EPS = 1e-12;

public:
    /*----------------array----------------*/
    static int FindEqualIndex(const std::vector<double>& vec, double value)
    {
        if (vec.empty() || IsLess(value, vec.front()) || IsGreater(value, vec.back())) return -1;
        int index = std::lower_bound(vec.begin(), vec.end(), value) - vec.begin();
        return IsEqual(value, vec[index]) ? index : -1;
    }

    static int FindNearestIndex(const std::vector<double>& vec, double value)
    {
        if (vec.empty()) return -1;
        if (value <= vec.front()) return 0;
        if (value >= vec.back()) return int(vec.size()) - 1;
        int idx = int(std::upper_bound(vec.begin(), vec.end(), value) - vec.begin()) - 1;
        double m = (vec[idx] + vec[idx + 1]) / 2;
        return value <= m ? idx : idx + 1;
    }

    /*----------------algebra----------------*/
    static bool IsValid(double x) { return std::isfinite(x); }

    static bool IsZero(double x, double eps = ZERO_EPS) { return x >= -eps && x <= eps; }

    static bool IsZero(const double* x, size_t n, double eps = ZERO_EPS)
    {
        if (!x) return true;
        for (size_t i = 0; i < n; ++i)
            if (!IsZero(x[i], eps)) return false;
        return true;
    }

    static bool IsZero(const std::complex<double>* x, size_t n, double eps = ZERO_EPS)
    {
        if (!x) return true;
        for (size_t i = 0; i < n; ++i)
        {
            if (!IsZero(x[i].real(), eps) || !IsZero(x[i].imag(), eps)) return false;
        }
        return true;
    }

    static bool IsZero(const int* x, size_t n)
    {
        if (!x) return true;
        for (size_t i = 0; i < n; ++i)
            if (x[i] != 0) return false;
        return true;
    }

    static bool IsEqual(double x, double y, double eps = ZERO_EPS) { return IsZero(x - y, eps); }
    static bool IsLess(double x, double y, double eps = ZERO_EPS) { return x < y - eps; }
    static bool IsGreater(double x, double y, double eps = ZERO_EPS) { return x > y + eps; }

    static bool IsVecLess(const std::vector<int>& x, const std::vector<int>& y)
    {
        assert(x.size() == y.size());
        for (size_t i = 0; i < x.size(); ++i)
            if (x[i] != y[i]) return x[i] < y[i];
        return false;
    }

    static bool IsIdentityTransMatrix(const double v[3][3])
    {
        for (int i = 0; i < 3; ++i)
            for (int j = 0; j < 3; ++j)
                if (i == j)
                {
                    if (!IsZero(v[i][j] - 1.0)) return false;
                }
                else
                {
                    if (!IsZero(v[i][j])) return false;
                }
        return true;
    }

    static double Norm2(const std::vector<double>& x) { return Norm2(x.data(), x.size()); }
    static double Norm2(const double* x, size_t n)
    {
        double r = 0.0;
        for (size_t i = 0; i < n; ++i) r += x[i] * x[i];
        return std::sqrt(r);
    }
    static double Norm2(const std::vector<std::complex<double>>& x) { return Norm2(x.data(), x.size()); }
    static double Norm2(const std::complex<double>* x, size_t n)
    {
        double r = 0.0;
        for (size_t i = 0; i < n; ++i) r += x[i].real() * x[i].real() + x[i].imag() * x[i].imag();
        return std::sqrt(r);
    }

    static double DotProduct(int n, const double* x, const double* y)
    {
        double r = 0.0;
        for (int i = 0; i < n; ++i) r += x[i] * y[i];
        return r;
    }
    static double DotProduct(int n, const std::complex<double>* x, const std::complex<double>* y)
    {
        double r = 0.0;
        for (int i = 0; i < n; ++i) r += std::abs(x[i] * y[i]);
        return r;
    }
    static double Dot(const double* x, const double* y) { return x[0] * y[0] + x[1] * y[1] + x[2] * y[2]; }
    static void Cross(const double* x, const double* y, double* r)
    {
        r[0] = x[1] * y[2] - x[2] * y[1];
        r[1] = x[2] * y[0] - x[0] * y[2];
        r[2] = x[0] * y[1] - x[1] * y[0];
    }

    static double Random()
    {
        std::srand(unsigned(std::time(nullptr)));
        std::rand();  // skip first
        double v = std::rand() / double(RAND_MAX);
        const double eps = 1e-4;
        if (v > 1 - eps) v = 1 - eps;
        if (v < eps) v = eps;
        return v;
    }

    static double Conj(double x) { return x; }
    static std::complex<double> Conj(std::complex<double> x) { return std::complex<double>(x.real(), -x.imag()); }

    /*----------------geometry----------------*/
    static double CalcDistance(const double* a, const double* b, int dim = 3)
    {
        std::vector<double> d(dim);
        for (int i = 0; i < dim; ++i) d[i] = b[i] - a[i];
        return Norm2(d);
    }
    static double CalcDistance2(const double* a, const double* b, int dim = 3)
    {
        double r = 0.0;
        for (int i = 0; i < dim; ++i)
        {
            double t = b[i] - a[i];
            r += t * t;
        }
        return r;
    }
    static double CalcPolygonArea2D(int n, const double* x, const double* y)
    {
        double r = x[0] * (y[1] - y[n - 1]);
        for (int i = 1; i < n - 1; ++i) r += x[i] * (y[i + 1] - y[i - 1]);
        r += x[n - 1] * (y[0] - y[n - 2]);
        return r * 0.5;
    }
    static double CalcTriangleArea(double x1, double y1, double x2, double y2, double x3, double y3)
    {
        return (x1 * (y2 - y3) + x2 * (y3 - y1) + x3 * (y1 - y2)) / 2.0;
    }
    static double CalcTriangleArea(double x[3], double y[3]) { return CalcTriangleArea(x[0], y[0], x[1], y[1], x[2], y[2]); }
    static double CalcQuadrangleArea(double x[4], double y[4])
    {
        return CalcTriangleArea(x[0], y[0], x[1], y[1], x[2], y[2]) + CalcTriangleArea(x[0], y[0], x[2], y[2], x[3], y[3]);
    }
    static double VectorDot(const double a[3], const double b[3]) { return a[0] * b[0] + a[1] * b[1] + a[2] * b[2]; }
    static bool Normalize(double v[3])
    {
        double len = Norm2(v, 3);
        if (len < 1e-60) return false;
        v[0] /= len;
        v[1] /= len;
        v[2] /= len;
        return true;
    }
    static bool PerpendicularProjection(const double vec1[3], const double vec2[3], double r[3])
    {
        double v[3];
        Copy(v, vec1, 3);
        if (!Normalize(v)) return false;
        double dot = VectorDot(v, vec2);
        for (int i = 0; i < 3; ++i) r[i] = vec2[i] - dot * v[i];
        Normalize(r);
        return true;
    }
    static double CalcSolidAngle(const double v1[3], const double v2[3], const double v3[3])
    {
        double m1 = Norm2(v1, 3), m2 = Norm2(v2, 3), m3 = Norm2(v3, 3);
        double d1 = VectorDot(v2, v3), d2 = VectorDot(v1, v3), d3 = VectorDot(v1, v2);
        double f1 = m1 * m2 * m3 + m1 * d1 + m2 * d2 + m3 * d3;
        double f2 = std::sqrt(2 * (m2 * m3 + d1) * (m1 * m3 + d2) * (m2 * m1 + d3));
        return f1 >= f2 ? 0.0 : 2 * std::acos(f1 / f2);
    }

    /*----------------allocate----------------*/
    static void Malloc(void*& p, size_t sz)
    {
        assert(p == nullptr);
        p = std::malloc(sz);
        if (!p) std::printf("Space not correctly allocated.");
    }

    static bool IsMultExceed(int a, int b)
    {
        int c = a * b;
        return (c / b) != a;
    }

    static std::vector<int> Seq(int count, int start = 0)
    {
        std::vector<int> v(count);
        std::iota(v.begin(), v.end(), start);
        return v;
    }

    /* ========================================================================*/
    template <typename T>
    static void Zero(T* x, size_t n)
    {
        std::memset(x, 0, n * sizeof(T));
    }

    template <typename T>
    static void Zero(std::vector<T>& x)
    {
        std::memset(x.data(), 0, x.size() * sizeof(T));
    }

    template <typename T>
    static void Copy(T* dst, const T* src, size_t n)
    {
        std::memcpy(dst, src, n * sizeof(T));
    }

    template <typename T>
    static void SafeCopy(T* dst, const T* src, size_t n)
    {
        if (n) std::memcpy(dst, src, n * sizeof(T));
    }

    template <typename T>
    static void Move(T* dst, const T* src, size_t n)
    {
        std::memmove(dst, src, n * sizeof(T));
    }

    template <typename T>
    static void Add(std::vector<T>& x, const std::vector<T>& y)
    {
        assert(x.size() == y.size());
        for (size_t i = 0; i < x.size(); ++i) x[i] += y[i];
    }

    template <typename T>
    static void Add(T* x, const T* y, long n)
    {
        for (long i = 0; i < n; ++i) x[i] += y[i];
    }

    template <typename T>
    static void Axpy(const T* x, T* y, T a, long n)
    {
        if (a == T(1))
            for (long i = 0; i < n; ++i) y[i] += x[i];
        else if (a == T(-1))
            for (long i = 0; i < n; ++i) y[i] -= x[i];
        else
            for (long i = 0; i < n; ++i) y[i] += a * x[i];
    }

    template <typename T>
    static void Axpy(const std::complex<T>* x, std::complex<T>* y, T a, long n)
    {
        if (a == T(1))
            for (long i = 0; i < n; ++i) y[i] += x[i];
        else if (a == T(-1))
            for (long i = 0; i < n; ++i) y[i] -= x[i];
        else
            for (long i = 0; i < n; ++i) y[i] += a * x[i];
    }

    template <typename T>
    static void Xpay(const T* x, T* y, T a, long n)
    {
        if (a == T(1))
            for (long i = 0; i < n; ++i) y[i] += x[i];
        else if (a == T(-1))
            for (long i = 0; i < n; ++i) y[i] = x[i] - y[i];
        else
            for (long i = 0; i < n; ++i) y[i] = x[i] + a * y[i];
    }

    template <typename T>
    static void Multiply(size_t n, T* x, T y)
    {
        if (y == T(1)) return;
        if (y == T(-1))
            for (size_t i = 0; i < n; ++i) x[i] = -x[i];
        else
            for (size_t i = 0; i < n; ++i) x[i] *= y;
    }

    template <typename T>
    static void Multiply(size_t n, std::complex<T>* x, T y)
    {
        if (y == T(1)) return;
        if (y == T(-1))
            for (size_t i = 0; i < n; ++i) x[i] = -x[i];
        else
            for (size_t i = 0; i < n; ++i) x[i] *= y;
    }

    template <typename T>
    static void Append(std::vector<T>& x, const std::vector<T>& y)
    {
        x.insert(x.end(), y.begin(), y.end());
    }

    template <typename T>
    static void Append(std::vector<T>& x, const T* y, size_t n)
    {
        size_t old = x.size();
        x.resize(old + n);
        Copy(&x[old], y, n);
    }

    template <typename T>
    static void Append(std::set<T>& x, const std::set<T>& y)
    {
        x.insert(y.begin(), y.end());
    }

    template <typename T>
    static void Append(std::unordered_set<T>& x, const std::unordered_set<T>& y)
    {
        x.insert(y.begin(), y.end());
    }

    template <typename ForwardIterator, typename T>
    static int FindIndex(ForwardIterator first, ForwardIterator last, const T& value)
    {
        return int(std::lower_bound(first, last, value) - first);
    }

    template <typename T>
    static int LowerBoundIndex(const std::vector<T>& x, const T& v)
    {
        return int(std::lower_bound(x.begin(), x.end(), v) - x.begin());
    }

    template <typename T>
    static int UpperBoundIndex(const std::vector<T>& x, const T& v)
    {
        return int(std::upper_bound(x.begin(), x.end(), v) - x.begin());
    }

    template <typename T>
    static void SortAndDedup(std::vector<T>& vec)
    {
        std::sort(vec.begin(), vec.end());
        auto it = std::unique(vec.begin(), vec.end());
        vec.resize(it - vec.begin());
        vec.shrink_to_fit();
    }

    template <typename T>
    static bool Contains(const std::vector<T>& v, const T& val)
    {
        return std::find(v.begin(), v.end(), val) != v.end();
    }

    template <typename T>
    static bool Contains(const std::set<T>& s, const T& val)
    {
        return s.find(val) != s.end();
    }

    template <typename T>
    static bool Contains(const std::unordered_set<T>& s, const T& val)
    {
        return s.find(val) != s.end();
    }

    template <typename Key, typename Val>
    static bool Contains(const std::unordered_map<Key, Val>& m, const Key& key)
    {
        return m.find(key) != m.end();
    }

    template <typename T>
    static bool Intersects(const std::unordered_set<T>& s1, const std::unordered_set<T>& s2)
    {
        if (s1.size() > s2.size())
            for (const auto& v : s2)
                if (Contains(s1, v))
                    return true;
                else
                    for (const auto& v : s1)
                        if (Contains(s2, v)) return true;
        return false;
    }

    /*----------------quick sort----------------*/
    template <typename K, typename V>
    static void QuickSort(K* keys, V* values, size_t start, size_t end)
    {
        if (start + 1 < end)
        {
            size_t p = QuickSortPartition(keys, values, start, end);
            QuickSort(keys, values, start, p);
            QuickSort(keys, values, p, end);
        }
    }

private:
    template <typename K, typename V>
    static size_t QuickSortPartition(K* keys, V* values, size_t start, size_t end)
    {
        K pivot = keys[(start + end - 1) / 2];
        size_t i = start, j = end - 1;
        while (true)
        {
            while (keys[i] < pivot) ++i;
            while (keys[j] > pivot) --j;
            if (i < j)
            {
                std::swap(keys[i], keys[j]);
                std::swap(values[i], values[j]);
                ++i;
                --j;
            }
            else if (i == j)
                return i + 1;
            else
                return i;
        }
    }

public:
    /*----------------allocate----------------*/
    template <typename T>
    static void Malloc(T*& p, size_t n)
    {
        assert(p == nullptr);
        p = (T*)std::malloc(n * sizeof(T));
        if (!p) throw std::bad_alloc();
    }

    template <typename T>
    static void Free(T*& p)
    {
        if (p)
        {
            std::free(p);
            p = nullptr;
        }
    }

    template <typename T>
    static void Delete(T*& p)
    {
        if (p)
        {
            delete p;
            p = nullptr;
        }
    }

    template <typename T, typename... Args>
    static void Delete(T*& first, Args&&... rest)
    {
        Delete(first);
        Delete(rest...);
    }

    template <typename T>
    static void Clear(std::vector<T>& vec)
    {
        for (auto& x : vec) delete x;
        vec.clear();
    }

    template <typename T>
    static void Clear(std::list<T>& lst)
    {
        for (auto& x : lst) delete x;
        lst.clear();
    }

    /*----------------class type----------------*/
    template <typename Base, typename Derived>
    static bool IsA(const Derived* ptr)
    {
        return dynamic_cast<const Base*>(ptr) != nullptr;
    }
};
